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One approach to the calculation ol cross sections for infrared-safe observables in high energy 
collisions at next-to-leading order is to perform all of the integrations, including the virtual loop 
integration, by Monte Carlo numerical integration. In a previous paper, two of us have shown 
how one can perform such a virtual loop integration numerically after first introducing a Feynman 
parameter representation. In this paper, we perform the integration directly, without introducing 
Feynman parameters, after suitably deforming the integration contour. Our example is the A*'- 
photon scattering amplitude with a massless electron loop. We report results for N = 6 and N = 8. 



I. INTRODUCTION 

The calculation of cross sections in the Standard Model 
and its extensions at next-to-leading order in perturba- 
tion theory inevitably involves computing virtual loop 
Feynman diagrams. In a process in which 2 partons scat- 
ter to produce n partons we integrate over the momenta 
{pi, . . . ,Pn} of the n final state partons. In a leading or- 
der calculation, for each choice of {pi, . . . ,Pn}, we mul- 
tiply the desired measurement function for that point by 
the squared tree level matrix element evaluated at that 
point. At next-to-leading order, we need also the one 
loop matrix element times the complex conjugate of the 
tree level matrix element, plus the complex conjugate of 
this product. The standard method for this kind of calcu- 
lation involves computing the one loop matrix element as 
a whole: for the chosen {pi, . . . ,Pn}, we compute the in- 
tegral that represents the loop graphs, with infrared and 
ultraviolet subtractions as necessary. The method for cal- 
culating this integral typically involves representing the 
integral in terms of "master integrals," whose values are 
known. There has recently been very significant progress 
in developing this method [HEllS]. 

There is another possibility. The loop graph is an in- 
tegral over a loop momentum I. One can write the in- 
tegration over {pi, . . . ,pn} and I as a single integration, 
so that for each choice of {pi, . . . ,Pn} in a Monte Carlo 
style integration, we also choose a momentum I. Then 
we multiply the measurement function and the complex 
conjugate tree amplitude by the integrand of the loop 
amplitude evaluated at {{pi, . . . ,Pn}, I)- The loop am- 
plitude contains singular factors l/{{l — Qi)'^ + ie). These 
singularities can be partly avoided by deforming the in- 
tegration contour into the complex / space. Thus one 
needs to specify what the deformed integration contour 



is to be. 

There are some possible variations on this method. In 
one variation, the integral over the energy /° is performed 
analytically by closing the integration contour, leaving a 
three dimensional integration over I. In another vari- 
ation, the loop integral is re-expressed using Feynman 
parameters x, so that we have an integration over either 
I and X or just x. Then we integrate over a complex con- 
tour in X. In any of these variations, the loop integral 
is evaluated by numerical Monte Carlo integration along 
with the integration over {pi, . . . ,p„}, so we may refer 
to this as the numerical Monte Carlo method. 

The numerical Monte Carlo method was implemented 
in Ref. [4 for three jet observables in electron-positron 
annihilation. Here the loop integral is expressed as a 
three dimensional integration over I. The variant of the 
numerical Monte Carlo method in which the loop inte- 
gral is expressed as an integration over Feynman param- 
eters X has been successfully implemented in NLO cal- 
culations by Lazopoulos, Melnikov and Petriello [5] for 
pp ZZZ -\- X and by Lazopoulos, McElmurry, Mel- 
nikov and Petriello [B] for pp ttZ -\- X. Furthermore, 
Anastasiou, Beerli and Daleo have used this method to 
compute the two loop amplitudes needed for gg ^ h me- 
diated by a heavy quark or a scalar quark [7]. In these 
methods the infrared singularities are eliminated by the 
method of sector decomposition . For the contour de- 
formation in X space, these authors follow the prescrip- 
tion given in Ref. [S]. 

In Ref. [S] , two of the present authors explored the nu- 
merical Monte Carlo method using the Feynman parame- 
ter representation by taking as an example the amplitude 
for 7 + 7 ^ (A^ — 2)7 through a (massless) electron loop, 
as depicted in Fig. [T] The integral to be evaluate has the 
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FIG. 1: Feynman diagram for the A''-photon amplitude. 
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where N{1) is the numerator function. This particular 
amphtude is useful as a test case because the loop inte- 
grand has infrared singularities, but (for a generic choice 
of {pi, . . . ,p„}) these singularities are integrable. Thus 
infrared subtractions are not needed and one can test 
the integration method without confronting a subtrac- 
tion method at the same time. The integral Eq. ^ was 
re-expressed using Feynman parameters x, giving an in- 
tegral with the form 
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where £ is a translated and Wick rotated loop momen- 
tum, ^ • ^ is the euclidian square of i, and A^(a:) is a 
quadratic function of the Feynman parameters x. The 
Feynman parameters are integrated over a certain de- 
formed contour C. 

Although the numerical Monte Carlo integral ^ is 
guaranteed to converge as the number of integration 
points becomes arbitrarily large, it needed to be proved 
that the method could produce the right answer with 
a practical number of integration points. For this rea- 
son Ai was evaluated at specified points {pi, . . . ,p„} by 
numerical Monte Carlo integration. The practicality test 
showed that one could obtain results for N = 6. For some 
helicity choices, the six gluon amplitude was known an- 
alytically [10 and the numerical results agreed with the 
analytical answer. For other helicity choices, the results 



have been confirmed by independent analytical calcula- 
tions by Binoth, Heinrich, Gehrmann and Mastrolia [lT| 
and by Ossola, Papadopoulos and Pittau [H]. 

In this paper, we ask whether one could use Monte 
Carlo numerical integration to perform the integration 
in Eq. ([T|) directly, without transforming the integral 
into the Feynman parameter form ([2]). In a practical 
application for which infrared subtractions are needed, 
as in Refs. [3 El E] , one would lose the possibility of con- 
structing the subtractions with the sector decomposition 
method. On the other hand, it is known |13j how to con- 
struct the subtractions directly in momentum space in 
a fashion that is analogous to that used for real emis- 
sion diagrams. If one were to perform the integration in 
Eq. (ij directly, there would be a certain advantage of 
simplicity. Additionally, one would avoid having a de- 
nominator A(a::) raised to a high power, which can create 
numerical convergence difficulties. 

The real apparent advantage of the Feynman parame- 
ter integration is that A(a;) is simply a quadratic function 
of X, so that it is easy to find a contour deformation that 
keeps us away from the zeros of the denominator. With 
Eq. ([ij , the denominator is a product of factors, each of 
which vanishes on a different surface in I space. For this 
reason, it is not immediately evident how to deform the 
integration contour in I. In the subsequent sections, we 
lay out a contour deformation with the required prop- 
erties. We then apply the same practicality test that 
was used for the form ([2|. We find that the direct form 
([ij gives us somewhat better numerical convergence than 
does the Feynman parameter form We do not know 
if the direct form is better in real applications. Indeed, 
it may well prove more useful to use the style of calcula- 
tion in which the integral M is evaluated as a whole in 
terms of master integrals. Alternatively, it may be that 
one method is better for some applications while other 
methods work better in other applications. We offer a 
method for performing the integral ([l]) directly in this 
paper in order to extend the range of available choices. 



II. INTEGRATING ON A DEFORMED 
CONTOUR 

The integrand in Eq. ([ij has singularities on the sur- 
faces {I — Qi)^ — 0. In order to avoid these singularities, 
we can deform the integration contour so that the loop 
momentum has an imaginary part. Call the complex loop 
momentum £ and let 
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Here and k'^ are the real and imaginary parts of l'^ 
and K is a function of I. With this notation, the integral 
is 
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In moving the integration contour we make use of the 
multidimensional version of the widely used one dimen- 
sional contour integration formula. A simple proof is 
given in the second paper of Ref. [i]. The essence of 
the theorem is that we can move the integration contour 
as long as we start in the direction indicated by the +iO 
prescription and do not encounter any singularities of the 
integrand along the way. 

In order to see what is required for the deformation, 
we consider a family of deformations, specified by 

«^(0 = A(0<(0 . (5) 

We imagine starting with an infinitesimal A and then 
increasing it to its final value, X[{1). Thus we consider 
< X{1) < Af(/). The denominator corresponding to 
propagator i is 

{l-Q, + = {l- Q,f - A(0' ^o{l)^ 

+ 2iA(0(Z-g.)-Ko(0 ■ 

Our first requirement is that we start the deformation 
in the direction specified by the -l-iO prescription. This 
means that on any of the of the surfaces (l — Qi)^ — we 
have {I — Qi) ■ ko(Z) > 0. 

Consider a a point I on the cone {I — Qi)^ — 0. The 
condition {I — Qi) ■ kq{1) > has a simple geometrical 
interpretation: that ko{1) points toward the interior of 
the cone. To see this, we consider the point I+Xkq, where 
A is infinitesimal and positive. Points on the interior of 
the cone have {I + Xkq — QiY > 0. Expanding to first 
order in A and using (l—Qi)^ = 0, we have 2XKQ-{l—Qi) > 
0. Similarly, the condition {I ~ Qi) ■ ko{1) = has the 
geometrical interpretation that ko(0 is tangent to the 
cone. 

We want to escape from the singularities if we can. 
This means that on {I — Qi)"^ = we would like {I — 
Qi) ■ «^o(0 > 0- This is easy as long as there is only one 
cone involved. We simply need to find a vector kq that 
points toward the interior of the cone. No deformation is 
possible for the point I = Qi. We cannot have (/ — Qi) ■ 
Ko(0 > if {l — Qi) = 0. Thus the point I = Qi is a pinch 
singularity, meaning that we cannot deform the contour 
to get away from it. 

At the intersection of two cones, deforming the contour 
is a little more subtle. Consider two cones {I — Qi)^ = 
and (l — Qj)^ = and suppose that K = Qj — Qi is a 
timelike vector. On the intersection of these two cones, 
we need a vector ko{1) that points towards the interior 
of both cones. It is geometrically evident that this is 
possible. If X is a spacelike vector, we also need a vector 
kq(1) that points towards the interior of both cones. It is 
also geometrically evident that this is possible. 

If K = Qj — Qi is a. lightlike vector, the cones meet 
along a line I — Qi — xK. If x > 1 or a; < 0, the inside 
of one of the cones is inside of the other, so that there is 
a range of vectors k,o{1) that point toward the interior of 
both cones. However if < x < 1, the inside of one cone 
is outside of the other, so that there is no vector kq{1) 



that points toward the interior of both cones. Specifically, 
for I = Qi + xK = Qj - (1 - x)K with < x < 1, we 
need xK ■ ko{1) > and (1 - x)K ■ kq{1) < 0. The best 
that we can do is have xK ■ kq{1) = (1 — x)K ■ kq{1) — 0. 
Since — 0, this is possible with kq{1) — c{x)K. Thus 
the contour is pinched: we can deform along K, but not 
in any other direction. With this deformation, we do 
not escape from the singularity. This pinch singularity is 
called the coUinear singularity. 



III. PREVIEW OF THE DEFORMATION 

In the subsequent sections, we define a contour de- 
formation quite precisely. In this section, we provide 
an informal statement of the main idea. Consider two 
cones {I — Qi)^ = and {I — Qj)^ = and suppose that 
K = Qj — Qi is a timelike vector with > 0. Let 

Ko = -cil-Q^) . (7) 

The coefficient c is a non-negative function of I . We want 

Ko-il~ Qi) > (8) 

on the surface {I — Qi)^ = 0. But 

I^0■{l-Q^) = -C{l-Q^)^ , (9) 

SO this requirement is automatically met. We also want 

Ko-{l~ Qj) > (10) 

on the surface {I — Qj)^ = 0. There are two cases to 
consider. For the backward light cone from Qj, we simply 
note that the point Qi lies inside this backward light cone. 
Thus for I on the backward light cone from Qj , kq points 
to the interior of this light cone and KQ-{l — Qj) > 0. This 
is illustrated in Fig. [2j For the forward light cone from 
Qj, Kq points in the wrong direction if c > 0. However, 
we can have kq ■ (l — Qj) — by taking c = on the 
forward light cone from Qj. 

This construction also works if K is lightlike with K° > 
0. The only difference is that now the point I = Qi lies 
on the backward light cone from Qj. This means that 
K.o-{l—Qj) > on the backward light cone from Qj except 
along the line from Qi to Qj, where kq ■ {I — Qj) = 0. 
That is, we escape from the singularity everywhere except 
along the coUinear pinch singular line. This is illustrated 
in Fig. d 

We see that the contour deformation of Eq. ([t]) achieves 
Ko • - Qj) > on (/ - Qj)^ = and kq ■ {I - Qi) > 
on {I ~Qi)^ — 0. What we need is to replace the > signs 
by > signs on the light cones except where the contour is 
pinched. We can achieve this by including several terms 
of the form suggested by Eq. Q plus, in certain regions 
of I, some terms pointing in a fixed timelike direction. 
We specify this choice in subsequent sections. 
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FIG. 2: Direction of the deformation kq = —c {I — Qi) for 
selected points on the backward hght cone from Qj when 
Qj — Qi is a timelike vector with a positive time component. 
The arrows, which represent the direction of kq, point to the 
interior of the backward light cone from Qj, so that kq • (/ — 
Qj) > 0. 




FIG. 3: Direction of the deformation kq = —c{l — Qi) for 
selected points on the backward light cone from Qj when Qj — 
Qi is a lightlike vector with a positive time component. For a 
generic point on the backward light cone from Qj, the arrows, 
which represent the direction of kq, point to the interior of 
the light cone, so that kq ■ {I — Qj) > 0. On the lightlike 
line along the intersection of the two cones, kq is parallel to 
{I — Qj), so that Ko ■ {I ~ Qj) = 0. 

IV. DEFINING THE INTEGRAL 

As in Ref. j9], we wish to compute the amplitude in 
quantum electrodynamics for scattering of two photons 
to produce N — 2 photons by means of a (massless) elec- 
tron loop. This process is illustrated in Fig. [l] The 
integral is given in Eq. ([T]). Electron line n in the loop 
carries momentum I — (5„, where Q„ is fixed and we in- 
tegrate over I. The momentum carried out of the graph 
by external photon n is 

P'li = Qn+l ^ Qn I (11) 

with = 0. 

We attack this simple problem because infrared 
counter terms are not needed: the original integral is 
infrared finite. The propagator denominators provide 
factors that would lead to logarithmic divergences after 
integration over the soft and coUinear regions. How- 
ever, these divergences are cancelled. For each elec- 
tron line there is a factor (/ — ^„). Thus the numer- 



ator provides a factor that removes the soft divergence 
from the integration region (/ — Q„) — > 0. Similarly at 
each vertex there is a factor (/- Q^+i) f-n^^n) (/- <^„), 
where e„(Pn) is the polarization vector of the photon. 
In the coUinear limit {I — Qn) — *■ xP, this gives a factor 
-x{l-x)fJ,^{Pr:)fn = -2x{l~x)fn e„(P„)-P„. This 
vanishes because e„(P„) • P„ = 0. Thus the numerator 
also provides a factor that removes each coUinear diver- 
gence. The loop integral is also finite in the ultraviolet 
as long as TV > 4. (For = 4 the integral is diver- 
gent by power counting, so an ultraviolet counter term is 
needed.) 

It will prove helpful to adopt a bit more notation. Two 
of the momenta Pi are the negatives of the momenta of 
the two incoming partons. We choose our labels so that 
these are P/v and Pa (for some A ^ N). We define 

P = -P/v and P = -Pa: 

P = Q.-Q.. ^^^^ 
P = Qa — Qa+i ■ 

We choose our reference frame so that the transverse mo- 
menta Pt and Pt of the two incoming particles vanish. 

There are some facts about the kinematics that can 
be easily be understood with the aid of Fig. |4] There 
we see the and 3 components of the momenta Qn for 
a sample event. The momenta Pi — Qi+i — Qi of the 
external particles join the points. In this illustration. Pi, 
^2, P3, and P4 represent final state particles. Then P5 is 
the negative of the momentum of an incoming particle. 
Next, Pg and P7 are the momenta of final state particles 
and Pg is the negative of the momentum of the other 
initial state particle. We can see from the figure without 
the need of an algebraic proof that Qb — Qa is inside the 
positive light cone if only outgoing particles are attached 
along the loop in the positive loop direction from line a 
to line b, as is the case, for instance for Q5 — Q2- We see 
also that Qi, — Qa is spacelike if precisely one initial state 
particle and at least one final state particle lie along the 
loop moving from a to b, as is the case, for instance, for 
Q5 — Qt- 

One can deform away from the singularities (l — Qi)^ — 
except where the contour is pinched along the straight 
lines I — Qi — xPi for < a; < 1. These are the lines 
depicted in Fig. |4] The endpoints of the line segments, 
where two line segments meet, are the soft singularities. 
The interiors of the lines, with < a; < 1, represent the 
coUinear singularities. At the soft singularities, I = Qi, 
the deformation has to vanish, 

no{Q,) = . (13) 

Along one of the coUinear lines, the deformation has to 
be parallel to the line, 

4{Q^+xP,)^c{x)Pl' , (14) 

for < x < 1. Thus we will have {l — Qi) ■ ko{1) = along 
one of the coUinear lines. Away from these lines we want 
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FIG. 4: Kinematics for the iV-photon amplitude, illustrated 
for N = 8. The sketch shows the and components of the 
loop momentum I. There are also two transverse components 
that come out of the plane of the paper and are not seen. 
The points are possible points I = Qi. The lines l — Qi— xPi, 
where Pi — Qi+i — Qi and < j: < 1, are also shown joining 
the points. The Pi are lightlike momenta. 



the denominators not to vanish. Thus we require 

{I - Qi) ■ Ko{l) > (15) 

on the surface {I — Qi)^ = except along the coUinear 
lines in this surface. We also demand that ko{1) be a con- 
tinuous function. Then Eq. ( 15 I also holds in a neighbor- 
hood of any point on {l — Qif^ excluding the coUinear 
lines. 



V. GEOMETRIC ARRANGEMENT OF THE 
LIGHT CONES 

Fig. [5]shows the kinematics as in Fig.[4]but with more 
information indicated. We show the P and components 
of the loop momentum I with the projections onto the l^- 

plane of the points I = Qn indicated. The projections 
of the light cones (I — QnY = ^i'^ the regions between 
the two dashed lines that pass through each Q„. Four 
shaded regions are indicated. 

Consider light cones with vertices in the left region, 
those with vertices {Qi, . . . , Qa]- The forward light cone 
from Qi is tangent to the backward light cone from Qi+i 
along the line from Qi to Qi+i- The forward light cone 
from Qi intersects the backward light cone from Qj for 
j > i + 1. The forward light cone from Qi is tangent 
to the forward light cone from Qi+i, while the forward 
light cones from Qj for j > z -I- 1 are nested inside the 
forward light cone from Qi. The backwards light cones 



FIG. 5: Kinematics for the A^-photon amplitude, illustrated 
for A'' = 8, showing light cones (/ — Qi)'^ — 0. In the illustra- 
tion, iV = 8 and A = 5. 



are similarly nested: the backward light cone from Qj is 
tangent to the backward light cone from Qj-i, while the 
backwards light cones from Qi for i < j ^ I are nested 
inside the backward light cone from Qj. 

The light cones with vertices in the right region, with 
vertices {Qa+i, • • ■ , Qn}i have analogous geometrical re- 
lationships with one another. 

The forward light cone from a vertex Qi on the left 
does not intersect a backward light cone from a vertex 
Qj on the right except for i — I, j = N, for which 
these cones are tangent along the line from Qi to Qn. 
The forward light cone from a Qi on the right does not 
intersect a backward light cone from a Qj on the left 
except for i = A + 1, j — A, for which these cones are 
tangent along the line from Q^ to Qa+i- 

The forward light cones with vertices Qi and Qn are 
tangent, as are forward light cones with vertices Qa+i 
and Qa- The other forward light cones with vertices 
in the left region intersect the forward light cones with 
vertices in the right region. These intersections are in the 
top region. 

The backward light cones with vertices Qi and Qn are 
tangent, as are backward light cones with vertices Qa+i 
and Qa- The other backward light cones with vertices 
in the left region intersect the backward light cones with 
vertices in the right region. These intersections are in the 
bottom region. 

We will make use of these geometrical properties in 
constructing the deformation. We note that it is possible 
to have A — 1 or A = N. In these cases, the picture looks 
rather different but the properties stated above still hold. 
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where 



FIG. 6; Illustration of the double parton scattering singular- 
ity. 



We can define coordinates x and x as follows: 
Qa+i) ■ P 



P P 
{l-Qi)-P 
P ■ P 



(16) 



Then the left region is a; < 0, a; > 0, the right region is 
X > 0, X < 0, the top region is x > 0, a; > 0, and the 
bottom region is a; < 0, i < 0. 



VI. DOUBLE PARTON SCATTERING 
SINGULARITY 

If the line from Qi to intersects the line from Qa+i 
to Qa, there is a pinch singularity at the intersection that 
gives a singular integral.^ This is the double parton scat- 
tering singularity. In general, these lines do not intersect, 
but they may be close to intersecting. 

Physically, the double parton scattering singularity 
arises when the process illustrated in Fig |6] is possible 
with on-shell intermediate states. Incoming photon N 
divides into a coUinear electron-positron pair. Incom- 
ing photon A divides into a coUinear electron-positron 
pair also. The electron from photon N collides with the 
positron from photon A so as to produce two or more 
outgoing photons, which must then have total transverse 
momentum equal to zero.^ The electron from photon A 
collides with the positron from photon N so as to pro- 
duce a different set of two or more outgoing photons, also 
with total transverse momentum equal to zero. Thus 
there is a singularity when the total transverse momen- 
tum of some subset of the final state photons, containing 
at least two photons, has total transverse momentum 
equal to zero. One is close to a double parton scattering 
singularity when <^ s. 

All points on the line from Qi to Qn have a; = 0. The 
point on this line that also has a; = is 



vi = ai Qi + (1 - ai) Q 



N 



(17) 



^ The amplitude as a function of the external momenta 
{pi, . . . ,Pjv} is non-analytic at this point. However, at least for 
N = 6, the amplitude does not become infinite as one approaches 
the singular point 1141 . 

^ Recall that we choose our reference frame so that the transverse 
momenta P and P of the two incoming photons vanish. 



O^x P- P 

= [aiQi + (1 - ai)QN - Qa+i] ■ P 
= [-aiP + QN-QA+i]-P , 



(18) 



so 



ai 

1 — ai = 



{Qn - Qa+i) ■ P 
P ■ P 

{Qa+i-Qi)-P 
PP 



(19) 



Similarly, all points on the line from Qa+i to Qa have 
a; = 0. The point on this line that also has a; = is 



vii = ail Qa+i + (1 - an) Q/ 



(20) 



where 



an 
1 - an 



(Q. 



P P 
(Qi - Qa+i) ■ P 



(21) 



PP 



The difference vi — vu is a spacelike vector. The separa- 
tion between the points is then measured by —{vi — vu)^. 
When vi — vu is small, the region in which the integrand 
is almost singular is near 



vi + Vu 



(22) 



In our computer code, we choose the origin of coordi- 
nates for / so that v — 0. However, our notation in this 
paper does not assume this choice. 



VII. THE DEFORMATION 

In this section, we propose a specific deformation func- 
tion kq(1). The final deformation k(Z) will be propor- 
tional to Ko(l), but with its size adjusted to ensure that 
it is not too large. We will first simply state the defini- 
tion of Ko(0- Then we will explain the rational for its 
various parts. 



A. The general formula 

The definition is 

N 

ko = -^c,(Z-Q,) + 5+(P + P)-5_(P + P) . (23) 
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The coefficients Cj and c± are non-negative functions of 
I. For the generic case 1<A<A^— Iwe define 

Cj = h^{l - Qj-i) h+{l - Qj+i) h^{l - Qn) 

X h+{l - Qa+i) 9(1) J e{2,...,A^l} , 
Cj ^ h^{l - Qj-i) h+{l - Qj+i) h^{l - Qa) 

xh+{l~Q,)g{l) J e{A + 2,...,N~l} , 
ci = h+{l ~ Q2) h_{l - Qn-i) h+{l ~ 5(0 7 

CA ^h^{l- Qa-i) h+{l - Qa+2) h-{l - Qn) g(0 , 
CA+i = h+{l - Qa+2) h-{l - Qa-i) h+{l - Qi) 9{l) , 
CM ^h^{l- Qn-i) h+{l - Q2) h-{l - Qa) 5(0 , 
c+ ^ h^{l - Qa) h^{l - Qm) 

X (x + x)d(x + x> 0) 9^(1) , 
c_ = h+{l - Qi) h+{l - Qa+i) 

X [-(a; + x)] e{x + x < 0) g+{l) . 

(24) 

We define the functions h±{l), 9±{l), and 9(1) below. For 
the special case A — 1, we define 

Cj ^ h^{l~ Qj-i) h+{l - Qj+i) h^{l - Qi) 

X /!+(;- Qi).g(0 J e {3,...,7V_l} , 
ci = h_{l - Qn-i) h+{l - Q3) 9{l) , 
C2 = h+{l-Q3)h+{l^Q,)g{l) , (25) 
CN ^ h^{l - QN-i)h^{l ~ Qi) g{l) , 
c+ = /i_(/ - Qn) (x + x) d{x + X > 0) 9-{l) , 
c- = h+{l - Q2) [-{x + x)] e{x + x< 0) 9+{l) . 

For the special case A — N — 1, we define 





= 




- Qj-i) h+{l ^ Qj+i) h^{l - Qn) 




X 


h. 


-{I^Qn)9{1) je{2,...,iV-2}, 


Cl 




-{I 


^Q2)h+{l~QN)9{l) , 






{I 


- QN-2)h-{l - Qn) 9{l) 7 


Cn 


= h. 


(l 


-Q2)h-{l~QN-2)9{l) , 


c+ 




{I 


- Qn-i) {x + x)e{x + x> 0) g-{l) 




= h. 


(l 


- Qi) [~{x + x)] e{x + x< 0) 5+(0 



(26) 



The various factors Cj contain factors h±{l — Qi) where 



h^{k) = 



\k\ + Eu)^ 



and 



h+{k) 



{\k\ + Eu)^ + Ml 
(\k\-Ek)^ + Aq 



^(Ek>-\k\) 



Ek < \k\ 



(27) 



(28) 



Here Ek and k are the energy and three-vector parts of 
fc in a frame in which P + P points along the time axis. 
Thus h-{k) = for /c e ^-(O) and h+{k) = for A: S 
C+(0). These functions depend on the parameter Mi 
with default value A'h = 0.05 [P-P]^/^. 



We include in cj a factor g{l), 
9{l) 
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{P - t;0)2 + (/ - t;)2 + M| 



(29) 



Here v is the vector defined in Eq. ( 22 1 , 71 is a dimen- 
sionless parameter with default value 71 = 0.7, and M2 
is a parameter with dimension of mass with default value 
M2 = [P-P]i/2. 

In c+ we include a factor 9-{l) and in c_ we include a 
factor 9+{l), where 



9±il) 



72 



1 + (liS/w)' 



(30) 



where 72 is a dimensionless parameter with default value 
72 = 1 and 



E = r 



„0 



{I - v)' + Ml 



1/2 



(31) 



Here M3 is a parameter with dimension of mass with 
default value M3 = [P-P]i/2. 

Now let us examine what this formula does. We con- 
sider first the generic case, 1 < A < N — 1, for which 
Eq. (24 1 applies. At the end of this section, we discuss 



the cases A ~ 1 and A ~ N — 1. 

B. Cone notation 

In order to make the subsequent analysis more com- 
pact, we adopt a notation for cones. We denote the for- 
ward light cone from Qi, 

(/-g,)' = 0, {l-Q,)-iP + P)>0 , (32) 

by C+{Qi). Let us also denote by C+{Qi) the cone 
C+{Qi) together with its interior: 

il-Q,)^>0, il-Q,)-iP + P)>0 . (33) 

Similarly, we denote the backward light cone from Qi, 

{l-Q,)^^0, (z_Q^) .(p + p) <o , (34) 

by C-{Qi) and this cone together with its interior by 

{l-Q^)^>0, {l-Q^)-{P + P)<0 . (35) 

C. The coefficients Cj for j € {2, . . . , A - 1} 

Let us examine Cj for jS{2,...,A — 1}. The term in 
kq proportional to Cj is 

«j = -Cj {I -Qj) ■ (36) 
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We need 

■ {I - Qi) > (37) 
for I G C±{Qi) for any i other than j. We have 

Kj-{l-Q,) = -Cjil-Q^f + CjiQj-Q^)-il~Q,) . (38) 
For / £ C±{Qi), this becomes 

•(?-Q,)=c,(Q, -O,0-(^-Q.:) • (39) 

Consider first the case j < i < A. Then Qj E C^{Qi). 
For / e C^{Qi), we then have kj ■ (l — Qi) > 0, as required 
(since we define Cj so that it is non-negative). For / e 
C+{Qi), we have the "wrong" sign, Hj-{l — Qi) < 0, so we 
need to define Cj so that it vanishes for / E C+{Qi). We 
note that C+{Qi) C C+{Qj+i). Thus we ensure that Cj 
vanishes for I E C+{Qi) by including a factor /i+(Z— Qj+i) 
in Cj. 

Consider next the case 1 < * < j- Then Qj E C+{Qi). 
For I E C+{Qi), we then have Kj-{l—Qi) > 0, as required. 
For I E C-{Qi), we have the "wrong" sign, kj ■ {I — Qi) < 
0, so we need to define Cj so that it vanishes for I E 
C_(Qi). We note that C^{Qi) C C^{Qj^i). Thus we 
ensure that Cj vanishes for I E C-{Qi) by including a 
factor h-{l — Qj-i) in Cj. 

Finally, consider the case A + 1 < i < N. In this 
case, Qj — Qi is spacelike, so that Kj ■ {I — Qi) can have 
either sign for I E C±{Qi). Thus we need to define Cj so 
that it vanishes for / E C±{Qi). We note that C+{Qi) C 
C+{Qa+i) and C-{Qi) C C+{Qn)- Thus we ensure that 
Cj vanishes for I E C±{Qi) by including a factors h^{l — 
Qa+i) and h^{l - Qn) in Cj. 

We also include a factor g{l) in Cj. The purpose of the 
factor g{l) is to turn off the deformation when I and /° 
are large. 

D. The coefficient Cj for j = A. 

We have discussed Cj for j E {2, . . . , A — 1}. The case 
j = A is similar, but not identical. 

There are no cones C±{Qi) for j < i < A. Therefore 
we do not need a factor h^{l — Qj+i) in Cj. We simply 
omit this factor. We do include a factor h^{l — Qj^i), as 
before. 

For A + 2 < i < N, Qj — Qi is spacelike, so that Kj-{1 — 
Qi) can have either sign for / E C±{Qi). Thus we need to 
define Cj so that it vanishes for I E C±{Qi). We note that 
C+{Q,) c C+{Qa+2) and C_(Q,) C C+{Qn)- Thus we 
ensure that Cj vanishes for I E C±(Qi) by including a 
factors h^(l — Qa+2) and h^il — Qn) in Cj. 

¥ov i = A+l, Qj - Qi is lightlike, Qj E C+{Qi). For 
I E C+{Qi),we then have Kj-{l—Qi) > 0, as required. For 
I E C^{Qi), we have the "wrong" sign, Kj ■ {I ~ Qi) < 0, 
so we need to define Cj so that it vanishes for I E C-{Qi). 
We note that C-{Qi) C C{Qn)- Thus the factor h-{l~ 
Qn) in Cj ensures Cj vanishes for I E C-{Qi). 

We include a factor g{l) in Cj as before. 



E. Tlie coefficient Cj for j — 1. 

The case j = 1 is essentially the same as the case 
j — A. We do not need a factor h^{l — Qj-i) in Cj 
because there are no cones C{Qi) for I < i < j- We 
ensure that Cj = for ^ e C±{Qi) forA+l<i<A^-l 
and for / E C+{Qn) by including factors h+{l — Qa+i) 
and h^(l — Qn-i) in Cj. 

F. Tlie coefficients Cj for j G {A + 1, . . . , N} 

We define the coefficients Cj for j E {A + 1, . . . , N} 
according to the same pattern that we used for j E 
{1,...,A}. 

G. The coefficient c+ 

The deformations Hj — —Cj (l—Qj) have been arranged 
so that (l — Qi) ■ Hj > for I E C±{Qi). We in fact want 
{l~Qi)-J2 '^j > ^ except on the lightlike lines along which 
the contour is pinched. A straightforward but tedious 
analysis shows that we have achieved that objective in 
the left region and in the right region in Fig. [s] However, 
we have {I — Qi) • ^ Kj = on the cones C±{Qi) in a 
large part of the top region. This is because Kj = in 
the intersection of C+(Q2) and C+{Qa+2)- 

This deficiency is easy to fix. In the top region, on 
one of the cones C+(Qi) we have k ■ {I — Qi) > for any 
timelike k with > 0. In particular, one could take 
K (x P + P. We therefore add a term 

5+ (P + P) (40) 

to kq. We include in the coefficient c+ factors h^{l — Qa) 
and h^{l — Qn) so that c+ = on all of the backward 
light cones C-.{Qi), for which P + P is in the "wrong" 
direction. 

We also include in c+ factors (x + x) 9{x + x > 0) 
and g-{l)- The purpose of these factors is to control 
the deformation in the region of large momenta. Near 
the forward light cones C+{Qi) we have E ^ u. In this 
region, g_(/) ~ 1, while (x+x) grows with I. This gives us 
a big deformation that can keep the contour well away 
from these light cones. However, for any fixed I, {x + 
x) g-(l) ^ as — > oo, thus turning the deformation 
off. 



H. The coefficient c_ 

In the bottom region, on one of the cones C-{Qi) we 
have K ■ {I ~ Qi) > for any timelike k with < 0. In 
particular, one could take k cx — (P + P). We therefore 
add a term 

-c_(P + P) (41) 
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to Ko- We include in the coefficient c_ factors h+{l — Qi) 
and h^{l — Qa+i) so that c_ = on all of the forward 
hght cones C+(Qj), for which —{P+P) is in the "wrong" 
direction. 

We also include in c_ factors ~{x + x) 9{x + x < 0) and 
g-\-{l) that serve to control the deformation in the region 
of large momenta. 



I. The special case A = 1 

The preceding discussion has concerned the generic 
case 1 < A < N — 1. When A = 1, the situation is 
a little different. 

The definition of Cj for j G {3, . . . ,N — 1} follows the 
logic of the generic case. 

For ci, we note that Ki oc —{l — Qi) points in the right 
direction on C-{Qn) and on C-^-{Q2) but in the wrong 
direction on C-{Qi) for 2 < i < — 1 and on C+{Qi) 
for 3 < I < A^. Wc arrange that Ci vanishes on the cones 
C_(Qi) for 2 < i < iV - 1 and on C+(Q,) for 3 < i < 
by including factors h-{l — Qn-i) and h+{l — Q3) in ci. 

For C2, we note that K2 c>c —(I — Q2) points in the 
right direction on C_(Qi) and on C-{Qi) for 3 < z < iV. 
However, K2 points in the wrong direction on C+(Qi) and 
on C+{Qi) for 3 < i < iV. We arrange that C2 vanishes 
on the cones C+((5i) and on C+(Qi) for 3 < i < iV by 
including factors h+{l — Q\) and — Qs) in C2. 

For Cat , the same logic leads us to include factors /i_ {I— 
Qi) and h-{l — Qn-i) in cn- 

For c+, we note that P+P points in the right direction 
on all of the cones C_|_((5i) but in the wrong direction on 
all of the cones C-{Qj). We arrange that c+ vanish on 
all of the cones C_ {Qj ) by including a factor h-{l — Qm) 
in c+. We also include in c+ factors {x + x) 6{x + a; > 0) 
and g-{l), as before. 

The construction of c_ foUlows the logic used for c+. 



J. The special case A = N — 1 



can be sure that we have avoided all singularities except 
those that cannot be avoided because they are pinch sin- 
gularities. On the other hand, we come very close to the 
singularities, which is not good for the numerical con- 
vergence of the integration. Thus we need to make A as 
large as we can. Then the integration contour is far from 
the product of the real axes and it is not so easy to see 
how to guarantee that in deforming away from one cone 
{I + iXno — QiY = we do not run into another cone 
[l + iXKQ-Qjf = Q. 

To ensure that the integrand does not have any singu- 
larities as A varies from zero to its chosen positive real 
value, we can take A to be the smallest of a number, A^ (Z), 
defined for each propagator, a general choice, Ao(Z), to be 
fixed later, and a fixed, I independent, constant Ac (the 
default is Ac = 1): 



= min[Ac, Ao(0, min{Ai(/)}] 



(43) 



In order to examine what Aj ought to be, we consider the 
ith denominator, 

A = {l-Qi + i\KQf = {l-Qif + 2iXKo■{l-Q^)-\'^Kl . 

(44) 

We note that the function Di vanishes at values of A given 

by 



A = 



i Ko - {I - Qi) 



(45) 



Consider what happens if Kq {I — Qi)^ > [kq ■ {I — Qi)]"^- 
Then one of the poles crosses the real A axis when, as we 
vary I, ■ {I — Qi) crosses zero. The absolute value of 
the pole position is 



lAI 



• (/ - Qi)? + [i4 (i - Qi? - («o • {I - Qi))^] 



'(l-Qi)' 



The construction of the deformation when A = N — 1 
follows the logic of the case A= 1. 



VIII. HOW FAR TO DEFORM 



We now take the deformation to be Z — > Z -|- iK, where 



Akq 



(42) 



and A is a scalar normalization factor. We must ensure 
that the integrand does not have any singularities as A 
varies from zero to its chosen positive real value. 

It is not immediately obvious how to achieve this end. 
If we take a very small value of A then we have effectively 
an infinitesimal deformation. We have arranged that our 
deformation for small A is in the right direction, so we 



(46) 

There is a dangerous region [kq ■ {I — Qi)? <C — 
Qi)^. In this dangerous region, one of the two poles is 
near the positive real A axis. If we increase A from past 
the real part of the pole position, we come very near 
the pole. If Ko • {I — Qi) = 0, we pass through the pole. 
Thus we need to keep the actual value of A smaller than 
the real part of the pole position when the parameters 
(kq • {I — Qi),HQ {I — Qi)^) are in the dangerous region. 

Consider also the case Kq {I — Qi)^ < [kq ■ {I — Qi)?- 
Then Di vanishes at values of A given by 



A 



lio- {I - Qi) 



±^[Ko-il-Qi)?-Kl (l-Qi)^^ . 



(47) 
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Both poles are on the imaginary A axis and there is no 
limitation on how big a positive real value value one can 
take for A. In fact, choosing a bigger value for A moves 
us further away from the poles. 

We now define A^. We first define A^ for the region 2[ko- 
{I — Qi)]^ < Kq{1 ~ Qi)^, which includes the dangerous 
region [kq ■ {I — Qi)]^ <^ Kq{1 — Qi)^- In this region, to 
make sure that A^ is small enough, we let it be half the 
absolute value of the pole position. That is 



A - 



for 2[ko • {I 



(48) 



For the region < Kq ~ Qi)^ < 2[ko • H - Qi)]'^, we 



define 



2 



Hq{1 Qi)'' 



for Q<kI{1- Q,f < 2[ko ■ {I - Q^) 



(49) 



Note that the definitions (48 1 and (49) match at Kq {I — 
Q^Y = 2[ko • {I - Qi)f. For the region „ Q^f < q, 
we define 



for kI {I - Q,Y < . 



(50) 



Note that the definitions (49 1 and (50) match at Kq (Z — 

Q^Y = 0■ 

This method needs some special consideration for the 
case of Xi when I approaches one of the coUinear singu- 
larities along the two lines I = Qi + z{Qi^i — Qi) and 
I = Qi + z(Qi_i — Qi) that meet at I = Qi. In these 
cases, the numerator in the definition of Xi vanishes in 
the limit and also the denominator, Kq, vanishes in the 
limit. Potentially, then, Xi approaches zero in the limit 
and its derivatives with respect to I are singular. This 
could lead to a singularity in the jacobian for the defor- 
mation, the determinant of the matrix 6^ + idn'^/dl'^. 

To analyze this, we write 



Kq 



where 



N 

- Qj) + c+{P + P)-~C-{P + P) 

C{l-Q,) + R, , 



(51) 



N 



N 



(52) 

We note that 

no-{l-Qi)^~C{l~Qi)^ + {l~Qi)-R, , (53) 



while 

Kl = C^{l-Q,f ~2C{l-Qi)-R, + R^, . (54) 
We can eliminate {I ~ Qi) ■ Ri between these to obtain 



(/ - Q.f = -l^nl-l,n^-{l-Q,) + i:, Hf 



(J2 



C 



(J2 



(55) 



So far, this is general. If we now consider the approach to 
the coUinear singularity, we have i?f 0. Furthermore, 
Ri approaches zero fast enough, like [(/ — QiYY, that we 
can set it to zero here. Then 



2ko ■ (/ - Q.) w ~C{1 - Q^Y 



1 

C 



(56) 



Using this relation, we can see that when I is close to 
the coUinear singularity, it is never in the region 2[kq ■ 
{I - Qi)Y < kI{1 - QiY for which the choice (48 1 for A 
applies. Indeed, squaring Eq. (56) gives 



1 



+ ^[a-Q.)T + ^[-g]' (57) 
>kI{i-q,Y ■ 

When I is near the coUinear singularity, it can be in 
either of the regions for which the choices for A^ given 
in Eqs. (49) and (50) apply, depending on the sign of 
Kq (/ — Qif- If we use Eq. (56) in Eq. (49), we have, for 
«o a - Qi? > 0, 



A = 



> 



nl{l-Q^Y + Cm-Q^??+HT/C' 



A{4Y 



(58) 



4C2 



If we use Eq. (56) in Eq. (50), we have, for k§ {I — QiY < 
0, 



A? 



c^[{i~Q.Y? + H]/c' 

4{nlY 



> 



1 

4C2 



(59) 



Thus near the either of the coUinear singularities along 
lines that meet ai I = Qi we have 



A ,: > 



1 

2C 



(60) 



We want to ensure that A is continuous as one ap- 
proaches the coUinear singularity. We define 



Ao(0 



1 



4C(0 



(61) 



Then we will have A = Ag = 1/(4 C), which is a smooth 
function of I, near the coUinear singularity. 
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IX. RESULTS 



We have constructed computer code [15, that carries 
out the calculations outlined here. In this section, we use 
this code to test how well the method described works. 
We calculate the sum over graphs oi A4 as specified in 
Eq. Q. For N photons, there are {N — 1)! graphs ob- 
tained by taking non-cycic permutations of the momenta 
Pi and photon polarization vectors e^.'^ We start with a 
labeling in which photons 1 and 2 are the incoming par- 
ticles and photons 3, . . . , iV are outgoing. Then we sum 
over graphs by summing over the results obtained with 
non-cyclic permutations of the indices i. 

The result for a given choice of helicities of the photons 
has a phase that depends on the precise definition of the 
photon polarization vectors e^. However, the absolute 
value of the scattering amplitude Ai is independent of 
the conventions used to define the e^, so we concentrate 
on Since \M \ is proportional to a^/^ and has mass 
dimension 4 — A^, we exhibit x {y/s)^~^ /a^^^ in our 
plots. We specify helicities in the form hi,h2,h^, ■ ■ ■ ,hiq , 
where 1 and 2 are the incoming particles and, following 
convention, hi and h2 are actually the negative of the 
physical helicities of the incoming photons. 

We compute M. by Monte Carlo numerical integration 
on the deformed contour as described in this paper. The 
integration samples points / giving extra weight to re- 
gions near the intersections of two or three light cones 
and to the region near the double parton scattering con- 
figuration. Of the {N — 1)! graphs, most weight is given 
to graphs for which the corresponding integral is clos- 
est to being singular. The code at [TF comes with some 
documentation of the sampling method. 

A. AT = 6 

We compute the six photon amplitude along a certain 
one dimensional curve in the space of final state mo- 
menta. We take photon 1 to have momentum pi along 
the — z-axis (so the physical incoming momentum is along 
the -I- z-axis), and we take p2 along the -t- z-axis. We 
choose an arbitrary point for the final state momenta 

{P3,P4,P5,P6}: 

P3 = (33.5,15.9,25.0) , 

P4 = (-12.5,15.3,0.3) , 

P5 = (-10.0,-18.0,-3.3) , ^ ' 

Pe = (-11.0,-13.2,-22.0) . 

Then we create new momentum configurations by rotat- 
ing the final state through angle 9 about the y-axis. 



With vector electromagnetic currents, charge conjugation invari- 
ance imphes that graphs that differ by reversing the order of the 
external photons are equal. Thus one might say that there are 
{N — l)!/2 independent graphs. 



30000 \ 




FIG. 7: Results for the six photon amplitude. An arbitrarily 
chosen final state was rotated about the y-axis through angle 

9. The numerical results for the helicity choice -I— I are 

compared with the analytic result of Ref. (TU]. The numer- 
ical results for the helicity choice H \—\ — are compared 

with the analytic result of Ref. pQ, . The numerical results 
were generated using 10^ Monte Carlo points for each of 120 
graphs. 



We compute the six photon amplitude by Monte Carlo 
numerical integration on the deformed contour as de- 
scribed in this paper. In Fig. [7] we plot computed values 
of s \ A4\/a^ versus 6 in the range from to tt. We show 

numerical results for the helicity choices +H and 

H hH — . For the helicity choice +-{ , we com- 
pare the results with the analytic result of Mahlon, |TU]. 
For the helicity choice H hH — , we compare the re- 
sults with the recent analytic result of Binoth, Heinrich, 
Gehrmann and Mastrolia ^llj, who also confirm the re- 
sult for +H . We see that the numerical results 

agrees with the analytical results. At « 2.32, the ex- 
ternal momentum configuration lies close to a double par- 
ton scattering singularity: {pt,3 +Pt,5)^ ~ 0.0003 s. We 
note that there is a quite sharp structure in \A4 \ near this 

angle for the helicity choice H 1 — I — . The numerical 

results were generated using 10^ Monte Carlo points for 
each of 120 graphs. 

Ref. shows analytically that Ai vanishes for the 
helicity choices ++++++ and ++++H — . We confirm 
that the numerical integration gives zero within errors 
for these helicity choices. 

Compared to the method of Ref. [S] that applies Monte 
Carlo numerical integration to the Feynman parameter 
representation of Eq. ([2]), how practical is the method 
presented in this paper? To see, we present in Fig. [8] 
numerical results using the Feynman parameter repre- 
sentation. The numerical results for < < 2.0 were 
published in Ref. [9j. For each of these angles, 1 x 10^ 
Monte Carlo points were used for each of 120 graphs. 
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FIG. 8: Results for the six photon amplitude using the Feyn- 
man parameter representation, Eq. (|2|. The labeling is as in 
Fig.0 



The points for 2.0 < < w were added later after the 
analytical results of Ref. [11] were published. For these 
angles, we used 3 x 10^ Monte Carlo points for each graph 
in order to explore with higher accuracy the region near 
6 = 2.3, which is close to a double parton scattering sin- 
gularity. The running time per Monte Carlo point in the 
Feynman parameter program is roughly half that of the 
direct momentum space method. Thus running times in 
Fig. [T] are comparable to those for Fig. [8] Comparing the 
error bars, we conclude that the direct momentum space 
method somewhat outperforms the Feynman parameter 
method. 



B. TV = 8 

We have also computed the eight photon amplitude 
along a one dimensional curve in the space of final state 
momenta, as for six photons. We take photon 1 to have 
momentum pi along the — z-axis and we take p2 along 
the -l-z-axis. We choose an arbitrary point for the final 
state momenta {p3,P4,P5,Pe,P7,P8}- 

P3 = (33.5,5.9,25.0) , 
Pi = (1.5,24.3,0.3) , 
P5 = (-19.1, -35.1, -3.3) , 
P6 = (28.2, -6.6, 8.2) , 
pV = (-12.2,-8.6,8.2) , 
PS = (-31.9,20.1,-38.4) . 



(63) 



Then we create new momentum configurations by rotat- 
ing the final state through angle 6 about the y-axis. 

We compute the six photon amplitude by Monte Carlo 
numerical integration on the deformed contour as de- 
scribed in this paper. In Fig. [9] we plot computed 



FIG. 9: Results for the eight photon amplitude. An ar- 
bitrarily chosen final state was rotated about the j/-axis 
through angle 6. The numerical results for the helicity 

choice -I— I are compared with the analytic result 

of Ref. |10| . For most of the chosen values of 0, the numerical 
results were generated using 2 x 10^ Monte Carlo points for 



each of 5040 graphs. For 6 = 2.0 and 9 ■ 
used for each graph. 



2.4, 10 points were 



values of s^|A^|/a^ versus 6 in the range from to 
TT. We show numerical results for the helicity choices 

H — I and compare the results with the analytic 

result of Mahlon [10] . There are 5040 graphs and we used 
2 X 10^ Monte Carlo points for each of them except at 
the angles 6 — 2.0 and 9 = 2.4, where we used 10^ points 
for each graph. We agree with the analytic results. As 
one can see from the figure, the numerical convergence 
of the integration is not as good with eight photons as 
it was for six. This should not be a surprise. There are 
now more graphs and, in each graph, there are now eight 
propagators that can be singular. The contour deforma- 
tion allows us to escape from singularities, but the more 
propagators there are, the less effective this escape is. 



Despite the marginal convergence of the integration 
with eight photons, we can note that the deformation 
method is quite robust. We have simply given it a hard 
test problem. The integration method of Ref. [9], in- 
volving numerical integration using Feynman parame- 
ters, could not produce any results at all with eight pho- 
tons. 
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X. CONCLUSIONS 



We have investigated how to perform the integral 



M 



dH 
J2ny 



N 



(-ie)^7V(0n 



(64) 



for the A^-photon scattering amphtude with a massless 
electron loop using numerical Monte Carlo integration. 
The main ingredient needed to perform the integration 
is to deform the integration contour into the space of 
complex momenta Z + in such a way as to avoid the 
singularities of the denominators wherever the contour is 
not pinched. 

We have in mind applications of the method described 
here to next-to-leading order calculations of infrared-safe 
observables in the Standard Model and its extensions. 



For this application, the integral of Eq. (64 1 is embed 



ded in a larger integral and the whole integral is per- 
formed by numerical Monte Carlo integration. Gener- 
ally, infrared subtractions are needed. These can be ob- 
tained for virtual loops from a construction jl3j that is 
similar to the construction commonly used for infrared 
singularities in real emission diagrams. The numerical 
integration method for virtual loop integrals has proved 
useful, as described in the Introduction, but it is not clear 
to us that this method is more practical than methods 
that calculate the amplitude Al as a whole and insert 
the complete value of M into the larger integral that fi- 
nally gives a physical observable. We offer this method 
because we believe that it is good to have a range of 
methods available. 

The method described here is similar to the method of 



Ref. [5] , in which the amplitude ( 64 1 is first represented as 
an integral over Feynman parameters. The Feynman pa- 
rameter method involves simpler contour deformations. 



On the other hand, our numerical results suggest that the 
convergence properties of the integral are better with the 
present "direct" method. 

We have seen how to deform the integration contour 
for the problem posed, with massless particles. This is an 
important case for applications to high energy scattering. 
It is also a case with some difficulties because with mass- 
less particles there are soft and coUinear singularities in 
the integrand. We leave for future work cases with par- 
ticle masses. With masses, the cones that appear in this 
paper turn into hyperboloids, requiring a different choice 
of contour deformations. Assuming that one has rede- 
fined the contour deformations, having non-zero masses 
can make the problem easier since it eliminates some or 
all of coUinear and soft singularities. On the other hand, 
the general problem with masses is more difficult because 
threshold singularities can occur. 

We also note that there could well be applications in 
which one wants to compute integrals like ( 64 ) as a whole. 



without needing to insert the integral into something else. 
The completely numerical method presented here has an 
advantage of simplicity: the method does not depend on 
the numerator function N{1). One might hope to have a 
general method that constructs the proper contour defor- 
mations for the denominators, with masses, as we have 
done here for the massless case. 
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